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USE OF THE TENSOR PRODUCT FOR NUMERICAL WEATHER 
PREDICTION BY THE FINITE ELEMENT METHOD - PART 1. 

Introduction 

In Ref. 1 Hinsman has developed a Finite Element (FE) 
program for Numerical Weather Prediction applications. The 
grid employed is rectangular with nodes at the intersections 
of north-south and east-west lines. It was shown by Stani- 
forth and Mitchell (Ref. 2) that the coefficient matrices 
for such a grid could be expressed as tensor products. In 
these products the factors are matrices which depend solely 
on grid spacing in the two orthogonal directions. This 
report deals with the coefficient matrix called the "mass" 
matrix in FE parlance. (In Refs. 3 and 4 applications of 
the tensor product resolution to the FE "stiffness" matrix 
are considered.) The theory which underlies the economical 
computational scheme based on the mass matrix resolution is 
first presented. Next, the number of floating point opera- 
tions and the number of storage locations needed for the 
coefficient matrix of this scheme are compared with those 
required by other better-known algorithms. A set of FORTRAN 
subroutines for implementing the tensor product scheme 
(TENSOR) is given in Appendix B. 



Consider the grid shown in Fig. 1. There are n spaces 



Theory 




10 11 12 Q 



2 ^3 ,4 in the E 



6 7 8 




rows 



(Note the cyclic 
boundary condition 
in the E-W direction.) 






X 

n columns 

Fig. 1. Node numbering and spacing. 



along each of e grid lines in the east-west direction. Node 
numbering is from west to east along successive grid lines, 

3 



beginning in the southwest corner. There is a cyclic bound- 
ary condition in the east-west direction so that the node 
number appearing at the beginning of each horizontal row is 
repeated at the end. Spacings of the horizontal and the 
vertical grid lines are not necessarily uniform. 



The computational problem addressed here is the solution 
of the equation 

Mw = V <1> 

where M (the "mass" matrix) is a square, symmetric matrix of 
size ne and w and v are column vectors of height ne. M and 
V are input quantities and w is sought. The tensor product 
representation of M is 

M = MB * MA <2> 

where MB is a square, symmetric, tridiagonal matrix of size 
e and MA is also square, symmetric and of size n. MA is 
tridiagonal except for nonzero elements in upper right and 
lower left corners. MB depends solely on the north-south 
node spacing b and MA depends upon the east -west node spac- 
ing a . The asterisk (*) denotes the tensor product. 
Explicit expressions for matrices MA, MB, and the tensor 
product are given in Appendix A. 

Let MB be represented as (e = 3) 



mb X X mb x 2 0 
MB = mbei mbaa mbas 
0 mb 3 2 mb 3 3 

If we partition w and v into e n x 1 subvectors so that 

n 

''I 

V * 



w = 



w 

w 



I 



II 

III 



III 



<4> 



we may use <2>, <3> and <4> to rewrite <1> as 

mbxi MA Wj + mbia MA Wjj = Vj 

mbai MA Wj + mb22 MA Wjj + mb 23 MA Wjjj - Vjj 

mbs 2 MA Wjj ^ mb 3 3 HA ~ '^III 
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Define 



H = <w, w 



I "ii "nP 



< 6 > 



and 




It is easy to verify that the equations <5> are equivalent 
to 



Solution of <7> can be accomplished by standard Gaussian 
elimination procedures. Specifically, the following steps 
are required. 



1) LDLT factoring of MA (n x n) . 

2) Forward reduction and back-substitution for e right- 



ilXCLllVX w ^ vJl w • 

(3) LDLT factoring of MB (e x e). 

(4) Forward reduction and back-substitution for n right- 
hand side vectors. 

This entire process is economical of both storage and arith- 
metic operations because of the tridiagonal structure of MA 
and MB . 

Boundary Conditions 

The cyclic boundary condition is implemented by repeating 
the node numbers of the western boundary on the eastern 
boundary as shown in Fig. 1. As already noted, this 
accounts for nonzero entries in the upper right-hand and 
lower left-hand corners of MA. 

It is sometimes required to impose a Dirichlet boundary 
condition on the southern and northern boundaries of the 
region. Specifically, the subvectors w and w of the solu- 
tion vector w are prescribed. To implement this boundary 
condition the following modifications to the standard 
solution procedure are required. 

In the n X e matrix V on the right-hand side of <7> the 
first and last columns are replaced by the prescribed bound- 
ary values of w, i.e., put Vj = Wj and Vg = w g. Let 
X = W MB and solve the system 



processing successive columns of V in standard fashion, but 
omitting the first and last columns. The reduced problem 



MA W MB = V 



<7> 



hand side vectors 



MA X = V 



< 8 > 
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now takes the form W MB = X and the first and last columns 
of X are Wj and Wg , respectively. Transposing both sides of 
this equation gives 

MB WT = XT <9> 

where WT and XT are the respective transposes of W and X 
(recall that MB is symmetric and is thus not altered on 
transposition) . Since the first and last rows of WT are 
known, the corresponding scalar equations are not needed. 
Accordingly, we form MBl by deleting the first and last rows 
of MB. We also reduce XT to XTl by omitting the first and 
last rows. This leaves the result 

MBl WT = XTl <10> 

or, in extenso, this takes the form (for n = 3, e = 5) 



mb 2 1 nib 2 2 nib 2 3 0 0 




k k iT 
u u u 




”k k k' 


0 mb 3 2 nib 3 3 mb 34 0 

0 0 mb43 mb44 mb 45 




u u u 
u u u 
_k k k_ 




k k k 
_k k k_ 



(In WT and XTl the elements denoted by "k" are known and 
those denoted by "u" are unknown. ) This equation may be put 
in standard form by first altering the first row of XTl by 
subtracting mb 2 i times the corresponding entries in the 
first row of WT and altering the last row of XTl by sub- 
tracting mb 45 times the corresponding entries in the last 
row of WT. Calling the new right hand side XT2 and forming 
MB2 from MBl by discarding the first and last columns and 
forming WTl by discarding the first and last rows of WT, the 
result is 

MB2 WTl = XT2 <11> 

Solution of <11> is carried out by LDLT factoring of MB2 , 
followed by forward reduction and back substitution. 



Floating Point Operations and Matrix Storage Requirements 

Presented here is a comparison of floating point opera- 
tion counts and matrix storage requirements for the tensor 
product scheme and three widely-used solution algorithms for 
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solving <7> or its equivalent <1>. Three of the schemes 
take advantage of symmetry of the coefficient matrices and 
store only elements on or above the principal diagonal. One 
of these, the **band solver" (BAND), places these elements in 
a rectangular matrix ne x r, where r is the maximum row 
length of the upper triangle of M. The "sky-line solver" 
(SKY) further economizes by storing only that part of the 
upper triangle beginning at the diagonal and extending to 
the topmost nonzero element of each column. These subvec- 
tors are assembled into a single vector. This scheme 
requires an additional integer address vector of length 
ne + 1. The remaining algorithm, "successive over-relax- 
ation" (SOR) , is iterative rather than direct. 

In most applications of the direct solvers the number of 
floating point operations required to factor the coefficient 
matrix into LDLT form is much greater than those required to 
complete the process of finding a single solution vector w 
corresponding to a given right-hand side vector v (forward 
reduction and back- substitution) . In the present applica- 
tion, however, the latter solution process must be carried 
out 17 times for each time step, so that the LDLT factoring 
makes a negligible contribution to the total computational 
expenditure. Accordingly, the operations required for fac- 
toring are not included in the tabulation below. 

In the following table the results given for the number 
of floating point operations are given in terms of the grid 
parameters n and e (defined in Fig. 1). One multiplication 
(or one division) plus one addition (or one subtraction) is 
counted as one operation. Exact results for these operation 
counts would take the form of a polynomial in n and e. Only 
the highest degree terms are given in the table. Since it 
is not possible to predict the number of iterations per 
solution when using SOR, the operation count given for that 
algorithm is for a single iteration . Also, since the number 
of storage locations required for SOR coefficient matrices 
is highly grid-dependent, no such entry is given for SOR. 
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TABLE I. Operation Counts and Storage Requirements. 



Algorithm 


Number of Operations 
per Solution 


Number of Storage Locations 
for Coefficient Matrices 


SOR 


10 en (1) 


(2) 


SKY 


2 en2 


en2 


BAND 


4 ep2 


2 en2 j 

1 


TENSOR 


8 en 


3 n + 4 e 



Notes: 1. Number of operations per iteration. 

2. Number of storage locations rs grid-dependent. 



Conclusion 

Close comparison of operation counts and storage require- 
ments leads to the conclusion that the TENSOR algorithm is 
clearly superior to the SKY and BAND algorithms. The com- 
parison with SOR is not as clear-cut. Considering, however, 
that the operation count for SOR is for only one iteration, 
there really seems to be little doubt that TENSOR is the 
method of choice. 
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APPENDIX A 

MATRICES MA, MB, AND THE TENSOR PRODUCT 



Symbols and b^. which appear in MA and MB are defined 
in Fig. 1. 



MA = I 

(n = 4) 



2(a4+Si) 0 

ai 2(ai+a2) 32 

0 32 2(32+83) 



34 



84 

0 

8 3 

33 2 ( 33 + 34 ) 



MB = f 
(e = 3) L 



2 bx bx 0 

bx 2(bx+b2) bx 
0 b 2 2 b 2 



The tensor product of matrices C and D may be represented 
in block partition form as 





"cxi D 


Cx2 0 


Cx 3 D 


C*D = 


C 2 X D 


C22 D 


C23 D 




.C 3 X D 


C32 D 


C33 0. 



where the c- • are the elements of C. Note that, if C and D 
* 

have dimensions r x s and t x u, respectively, the tensor 
product has dimensions rt x su . 
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APPENDIX B 

FORTRAN PROGRAM LISTINGS 



FORTRAN programs for the implementation of TENSOR are 
listed here. They appear in the form of subroutines AMTRX2, 
FACTOR, BACKA, and BACKB within the test program GAUSS3. 
(The subroutines FACTOR, BACKA, and BACKB are adapted from 
subroutine COLSOL of Ref. 5.) Also included are G0G3 , an 
Exec used to execute GAUSS2 - a dimensioned version of 
GAUSS3, and CDIM, an Xedit program used to enter the dimen- 
sions . 
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Listing: GAUSS3 FORTRAN 

MAIN PROGRAM MASS MATRIX USING TENSOR PRODUCT FACTORS 



THIS PROGRAM IS DESIGNED TO TEST THE SCHEME (TENSOR) 
WHICH RESOLVES THE MASS MATRIX INTO A TENSOR PRODUCT IN 
ORDER TO SOLVE THE SYSTEM OF EQUATIONS M w = v. THE 
SUBROUTINES MAY BE INSERTED IN THE PROGRAM DEVISED BY 
HINSMAN. 

:Z) 



IMPLICIT REAL*8CA-H.O- 
COMMON/ CMIA/NLAT , NLONC 

COMMON/ CMS /A(Z1),bTz1) , , . . . 

COMMON AG(ZB) ,BG(ZC) ,GA(ZK) ,GB1(ZL) ,GB2 (ZL) ,MA(ZM) , 
IMB(ZN) 

DIMENSION V(ZP) 

READ(5, *)NLONG, NLAT 
LATX=NIAT+1 
WRITEX6 ,1000) 

1000 FORMAT(//, MASS MATRIX - TENSOR PRODUCT RESOLUTION' 



^ItE (6,1001) NLONG , NLAT 
READ75,‘’')A.B 
WRITE(6 ,500)A 



503 
500 
1001 
C 



WRITEi6,5Q3)B 
FORMATU , B: 

FORMATf/. A: 

FORMAT ( ' NLONG 






NLAT = ' ,13, / ) 



501 

504 

1002 

1004 

1005 



WRITE! 
FORMA1 
WRITE! 
FORMA! 
WRITE! 
FORMA! 
WRITE! 
FORMA! 
WRITE I 
WRITE! 



,(12F4.1)) 
,(12F4.1)) 
(3X,6F7.3) ) 



! , (24F3. , 

..... - 

C0NSTRU"CT^ FACTORS V”gA, 6b 1’, AND GB2,"’OF I^S^'l^TRIX 
CALL AMTRX2 
WRITE (6, 501) AG 

^(i002)GA| 

^,i004)GBl’^ 

^ (603)m^^^ 

„*v^..v 6 ’l 006 )MB 

CU=GB1(3) 

CL=GB1X2*"'LATX-1) 

K=(NLAT-l)*NLONG 

C IF NDIR>0, THERE IS A DIRICHLET BOUNDARY CONDITION ON 
C NORTHERN AND SOUTHERN BOUNDARIES. 

READ(5,*)NDIR,V 

502 ]TORMaI(K^^1dIR^=^' , 11, / ' V: ' ,6F8.2,/ , (4X,6F8.2)) 
C PERFORM LDLT FACTORING OF GA, GBl, AND GB2 
CALL FACTOR (GA.MA, NLONG) 

CALL FACTOR (GBl, MB, LATX) 

CALL FACTOR (GB2, MB, LATX) 

WRITE (6, 1002 )GA 
WRITE (6, 1004 )GB1 
^6!i005)GB2 

ORWARD REDUCTION AND BACK- SUBSTITUTION USING 



C 

C 

C 

C 



C 

C 



C 

3 



C 

C 



ON NORTH AND SOUTH 



WRITE! 

PERFORM I 
FACTORS OF GA, 

CALL BACKA(GA,V,MA,NDIR) 

DIRICHLET BOUNDARY CONDITION 
BOUNDARIES ? 

IF(NDIR.GT.O)GO TO 3 
WRITE? 6, 5 10 W 

PERFORM FORWARD REDUCTION AND BACK- SUBSTITUTION USING 
FACTORS OF GBl 

CALL BACKB(GB1,V,MB,NDIR) 

GO TO 6 

CORRECT RIGHT-HAND SIDE FOR DIRICHLET CONDITION 
DO 2 J=l, NLONG 



V ( J + NLONd j = V ( J + NLONG ) - CU*Y (Jf } 

V( J + K ) = V 7 J + K ) - CL*V ( J + NLAT - NLO 



__ . __ NG) 

WftlTEf6,5lO)Y 

PERFORM FORWARD REDUCTION AND BACK -SUBSTITUTION USING 
FACTORS OF GB2 

CALL BACKB?GB2,V,MB,NDIR) 
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6 WRITE (6, 5 10 )V . 

510 FORMAT(/,’ V: ',6F8.2,/, 

C READ A NEW RIGHT-HAND SIDE 
READ(5,*)NDIR,V 
WRITET6,502)NDIR,V 
CALL BACKA(GA,V,MA,NDIR) 

if7ndir.gt.o)go to 7 

WRITE? 6, 5 low 
CALL BACKB(GB1,V,MB,NDIR) 
GO TO 16 
DO 5 J=l,NLONG 



nd’perfo^ a second solution 



7 

5 



16 

1003 

1006 



C ***• 

C . . 

c 
c 
c 
c 
c 
c 
c 
c 
c 



V ( J+NLONd )= V ( J+NLONG ) - CU*V U) 

vJj+K)=v7j+KJ- CL*V ( J+ NLAT'-NLONG ) 

WRITE? 6, 5 low 

CALL BACKB?GB2,V,MB,NDIR) 

WRITET6.516)V 

FORMAT?/, MA: ’ ,2X,36I3) 

FORMAT(/ ' MB: ’ ,2X,36I3) 

STOP 












SUBROUTINE FACTOR ( A, MAXA,NN) 



A(NWK) 

ma3La(n 



INPUT VARIABLES - - 

STIFFNESS MATRIX STORED IN COMPACTED FORM 
(MP) = VECTOR CONTAINING ADDRESSES OF DIAGONAL 
ELEMENTS OF STIFFNESS MATRIX IN A 
= NUMBER OF EQUATIONS 

NUMBER OF ELEMENTS BELOW SKYLINE 



NN 

NWK 

- - OUTPUT 
A (NWK) 



= D AND L - FACTORS OF STIFFNESS MATRIX 



C 

C 

C 

40 



50 



60 



70 

80 

90 



100 

110 

120 

140 

2000 



IMPLICIT REAL*8(A-H,0-Z) 

DIMENSION A(1),MAXA(1) 

PERFORM L*D*LT FACTORIZATION OF STIFFNESS MATRIX 

DO 140 N=1,NN 
KN=MAXA(N) 

KL=KN+1 

KU=MAXA(N+1)-1 

KH=KU-KL 

IF^KH^110,90,50 

IC=0 

KLT=KU 

DO 80 J=1,KH 
IC=IC+1 
KLT=KLT-1 
KI=MAXA(K) 

ND=MAXA(K+1)-KI-1 

IF(ND)80,80,60 

KK=MINO(ic,ND) 

DO 70 L=1,KK 

c=c+a(ki+l)*a(klt+l) 

A(KLT)=A(KLT)-C 

K=K+1 

K=N 

B=0. 

DO 100 KK=KL,KU 

K=K-1 

KI=I^XA(K) 

C=A(KK)/AfKI) 

B=B+C*A(KK) 

A(KK)=C 

AiKNl=A(KN)-B 

IF ( AiKN]) 120.120 , 140 

WRITE(I0UT,2000)N,A(KN) 

STOP 

CONTINUE 

FORMAT?/ / , ' STOP - STIFFNESS MATRIX NOT POSITIVE 
IDEFINITE. , // , ' NONPOSITIVE PIVOT FOR EQUATION 

END 
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C ** 






c 

c 

c 

c 



c 

c 

c 



c 

c 



140 

150 



160 



170 

180 

c 

c 

c 



SUBROUTINE BACKA(A,Y,MAXA,NDIR) 

4b*« mf^ 

^ ^ 4% ^ ^ r% ^ ^ r% ^ ^ #v ^ ^ 

THIS SUBROUTINE PERFORMS THE FORWARD REDUCTION AND BACK- 
SUBSTITUTION USING THE FACTORS OF GA 

IMPLICIT REAL*8(A-H,0-Z) 

COMMON/ CM lA/NLAT .NLONG 
DIMENSION A(l) ,V(1) ,MAXA(1) 

REDUCE RIGHT-HAND- SIDE LOAD VECTOR 

JMIN=1 

LATX=NLAT+1 

JMAX=LATX 

IS THERE A DIRICHLET BOUNDARY CONDITION? 

ifTndir.lt.iTgo to 140 
SKIP NORTH AND SOUTH BOUNDARIES 
JMIN=2 
JMAX=NLAT 

DO 240 J=JMIN, JMAX 
DO 180 N=l,NLONG 
KL=MAXA(N)+1 
KU=MAXA(N+1T-1 
IF^KU-KL) 180 , 160 , 160 

C=0. 

DO 170 KK=KL,KU 
K=K-1 

C = C + A ( KK ) *V ( K+ ( J - 1 ) *NLONG ) 

V j[N+ ( J - 1 ) *NL0NG> V ( N+ ( J - 1 ) *NLONG ) - C 
CONTINUE 



200 

210 

220 

230 

240 



C 

C 

C 

C 

C 

C 



C 

C 

C 



BACK- SUBSTITUTE 

DO 200 N=l,NLONG 
K=MAXA(N) 

V(N+(J-1) *NLONG ) = V ( N + ( J - 1 ) *NLONG ) / A ( K ) 

d 5 230 L=2,NLONG 
KL=MAXA7N)+i 

ku=maxa(n+iT-i 

IF^KU-KL)230 ,210,210 

DO 220 KK=KL,KU 
K=K- 1 

V(K+{J-l)*NLONG)=V(K+(J-l)*NLONG)-A(KK)*V(N+(J-l)* 
INLONG ) 

N=N-1 

CONTINUE 

RETURN 

END 









w ^ ^ ^ ^ ^ ^ ^ ^ ^ ^ ^ ^ ^ ^ ^ ^ ^ ^ ^ ^ ^ ^ 

. . . . S UBROUTINE. BACKB( A. V.MAKA.NDIR) 

THIS SUBROUTINE PERFORMS THE FORWARD REDUCTION AND BACK- 
SUBSTITUTION USING THE FACTORS OF GBl OR GB2 



IMPLICIT REAL*8(A-H.O-Z) 
COMMON/CMIA/NLAT .NLONG 
DIMENSION A(l) ,V(l) ,MAXA(1) 

REDUCE RIGHT-HAND- SIDE LOAD VECTOR 



LATX=NLAT+1 
NMIN=1 
NMAX = LATX 

C IS THERE A DIRICHLET BOUNDARY CONDITION? 

IFTnDIR.LT. l)GO TO 50 
C SKIP NORTH AND SOUTH BOUNDARIES 
NMIN=2 
NMAX=NLAT 

50 DO 240 J=l,NLONG 

150 DO 180 N=NMIN,NMAX 

KL=MAXA(N)+1 
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160 



170 

180 

C 

C 

C 



200 



205 



210 



220 

230 

240 



C 

C *• 



KU=MAXA(N+1)-1 
IF^KU-KL) 180 , 160 , 160 

C=0. 

DO 170 KK=KL,KU 
K=K-1 

C= C + A (KK ) *V ( J + f K - 1 ) *NLONG ) 
Vf J+?N-1)*NL0NG)=V(J+ ' ■ 

CONTINUE 

BACK- SUBSTITUTE 



(N-lj*NLONG)-C 



DO 200 N=NMIN,NMAX 

V?J+(N-1)*NL0NG)=V(J+(N-1)*NL0NG)/A(K) 

LMAX=LATX 

IFfNDIR.LT. 1)G0 TO 205 

LMIN=3 

LMAX=NLAT 

N=LMAX 

DO 230 L=LMIN,LMAX 

kl=maxa7n)+i 

KU=MAXA(N+1)-1 

IF(KU-KL)230,210,210 

d5 220 KK=KL,KU 
K=K- 1 

V(J+(K-1)*NL0NG)=V(J+(K-1)*NL0NG)-A(KK)*V(J+(N-1)* 

INLONG) 

N=N-1 

CONTINUE 

RETURN 

END 



C 

C 

C 

C 

C 

C 

C 

C 

C 

C 

C 

C 

C 



C 

C 



W 



2 

C 



SUBROUTINE AJMTRX2 

THIS SUBROUTINE FORMS THE MASS MATRIX IN THE FORM OF A 
TENSOR PRODUCT OF THE GB MATRIX AND THE GA MATRIX. 

THE FIRST OF THESE IS NLAT + 1 BY NLAT + 1, SYMMETRIC, 

AND TRIDIAGONAL. THE SECOND IS NLONG BY NLONG, 
SYMMETRIC, AND TRIDIAGONAL EXCEPT FOR SINGLE ELEMENTS 
IN UPPER RIGHT HAND AND LOWER LEFT HAND CORNERS. GB IS 
STORED IN SKYLINE VECTOR FORM TUPPER TRIANGLE WITH SPACE 
FOR FILL-IN) AS GBl AND GB2 . THE LATTER VERSION IMPOSES 
A DIRICHLET BOUNDARY CONDITION ON THE NORTH AND SOUTH 
BOUNDARIES. GA IS ALSO STORED IN SKYLINE VECTOR FORM. 
INTEGER ADDRESS VECTORS MB AND MA ARE ALSO GENERATED. 

IMPLICIT REAL*8 ( A-H . 0- Z ) 

COMMON / CMIA / NLAT , NLONG 

I^ESon^AGcIbS ,b1(z?)1’Ia(ZK) ,GB1(ZL) ,GB2(ZL) ,MA(ZM) , 

DIMENSION BG (NLAT) , AG (NLONG) ,GB1(2*NLAT- 1) , 

GA(3*NLONG-31 ,MA(NLONG+l) ,MB(NLAT+2) 

LATX=NLAT+1 

FIND BG = (ELEMENT HEIGHT) /6. 

DO 2 J=1,NLAT 

BG( J) =Bri+NLONG*( J- 1) ) / 3 . 

GENERATE GBl AND GB2 
GB1(1)=2.*BG(1) 

GB2(1)=1. 

DO 4 J=2,NLAT 
K=2"'( J-l) 

GB1(K)=2.*(BG(J-1)+BG(J)) 

GBlfK+l)=BG(jU) 

GB2CK)=GB1Tk) 

GB2(K+lJ=GBlfK+l) 

GB 1 ( 2*NLAT ) = 2 . ‘-BG (_NLAT ) 

GB 1 ( 2*NLAT + 1 ) = BG ( NLAT ) 

GB2(3)=0. 

GB2(2‘%LAT) = 1. 

GB212*NLAT+1)=0. 

FIND AG = (ELEMENT WIDTH) /6. 

15 
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DO 10 J=l,NLONG 

ag(j1=a(j)/3 . 

GENEME GA 

GA(l1=2.*(AG(l)+AG(NLONG)) 
DO 12 J=2>L0NG 
K=2*(J-1) 

GATK)=2.ffAG(J-l)+AG(J)) 
GA{K+ 1) — AG ( J ” 1 ) 

Kl=2*NLONG 
K2=3*NLONG-4 
DO 14 K=K1,K2 
GA(K)=0. 

GA] 3 '%L0NG - 3 ) = AG (NLONG ) 
GENERATE DIRECTORY VECTORS 
MB(1)=1 
D016 J=1.NLAT 
MBfj+l)=2*J 
MB C NLAT+ 2 ) = 2* ( NLAT+ 1 ) 

maU}=i 

DO 18 J=2,NLONG 



MA(j5=2*(i-l). 

MAiNLONG+l)=3* 



RETURN 
END 



NLONG-2 
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LiXStxng* GOGS EX£C 

ERASE GAUSS 2 * Al 

COPY GAUSS3 FORTRAN Al GAUSS2 = 

&STACK CDIM 

& STACK FILE 

X GAUSS 2 FORTRAN 

FORTGI GAUSS 2 

GLOBAL TXTLIB FORTMOD2 MOD2EEH 
FILEDEF 05 DISK 
LOAD GAUSS2 (START 



Listing: CDIM XEDIT 

SET CMSTYPE HT 
TOP 

* ZB=NLONG 
C^^ZB/6/ 

* ZC=NLAT, 

C^^ZC/3/ « 

* zi=nlong*ni;at 

5o^Z1/18/ ^ 

* ZK=3*NLONG-3 
C /ZK/15/ * * 

TOP 

ZL=2*NLAT+1 

U^W7/ 

* ZM=NLONG+l 
C^^ZM/7/ - - 

* ZN=NLAT+2 
yZN/5/ 

* ZP=NLONG*(NLAT+l) 

cJzP/24/ * * 

SET CMSTYPE RT 
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